clear;

% Switching in markup process

x=zeros(16,1);

x(1) = 2.0995;        % alpha1
x(2) = .96469;      % alpha2
x(3) = .08014;         % gamma1
x(4) = .11189;         % gamma2
x(5) = .997977;       % beta1
x(6) = .39149;        % kappa1;
x(7) =  3.6836;         % tau1;
x(8) = .09758;         % rho_a;
x(9) = .5;         % rho_f1;
x(10) = .9539;         % rho_f2
x(11) = .6656;          % rho_i
x(12) = .005;      % sd_a
x(13) = .005;      % sd_f1
x(14) = .0034;     % sd_i
x(15) = 0.00811;   % piess
x(16) = 0.006;     % g_y
x(17) = .00033;    % sigu
x(18) = 9.52;      % lnA_0
x(19) =-0.07;      % ystar
x(20) = .955377;       % p11
x(21) = .959815;      % p22
x(22) = .925;       % p11_f
x(23) = .9;       % p22_f

para = x;

alpha1 = para(1);
alpha2 = para(2);
gamma1 = para(3);
gamma2 = para(4);
beta1  = para(5);
kappa1 = para(6);
tau1   = para(7);
rho_a = para(8);
rho_f1 = para(9);
rho_f2  = para(10);
rho_i  = para(11);
sd_a  = para(12);
sd_f1  = para(13);
sd_i   = para(14);
piess  = para(15);
gy     = para(16);
sigu   = para(17);
lnA_0  = para(18);
yst    = para(19);
p11 = para(20);
p22 = para(21);
p11_f = para(22);
p22_f = para(23);

P_i = [p11 1-p11;
       1-p22 p22];
   
P_f = [p11_f 1-p11_f;
       1-p22_f p22_f];   

P = kron(P_i,P_f);

iss    = gy-log(beta1) + piess;

%----------------------------------------------------------------------
% check the roots of the system using gensys
%----------------------------------------------------------------------
A = [p11 1-p11 (p11/tau1) ((1-p11)/tau1);
     1-p22 p22 ((1-p22)/tau1) (p22/tau1);
     0 0 beta1*p11 beta1*(1-p11);
     0 0 beta1*(1-p22) beta1*p22];

B = [1+(gamma1/tau1) 0 (alpha1/tau1) 0;
     0 1+(gamma2/tau1) 0 (alpha2/tau1);
     -kappa1 0 1 0;
     0 -kappa1 0 1];


if (min(abs(real(eig(inv(A)*B))))) < 1
    disp('*************  DL determinacy condition violated  *************')
end


%----------------------------------------------------------------------
% compute the solutions for inflation and output
%----------------------------------------------------------------------
A_Aa_mat = [1-p11*rho_a, -(1-p11)*rho_a, -(p11*rho_a)/tau1, -((1-p11)*rho_a)/tau1, 1/tau1, 0;
            -(1-p22)*rho_a, 1-p22*rho_a, -((1-p22)*rho_a)/tau1, -(p22*rho_a)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_a, -beta1*(1-p11)*rho_a, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_a, 1-beta1*p22*rho_a, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];           
        
solvec_a = inv(A_Aa_mat) * [rho_a/tau1;rho_a/tau1;0;0;0;0];


A_Af_mat_11 = [1-P(1,1)*rho_f1, -P(1,2)*rho_f2, -P(1,3)*rho_f1, -P(1,4)*rho_f2;
               -P(2,1)*rho_f1, 1-P(2,2)*rho_f2, -P(2,3)*rho_f1, -P(2,4)*rho_f2;
               -P(3,1)*rho_f1, -P(3,2)*rho_f2, 1-P(3,3)*rho_f1, -P(3,4)*rho_f2;
               -P(4,1)*rho_f1, -P(4,2)*rho_f2, -P(4,3)*rho_f1, 1-P(4,4)*rho_f2];
           
A_Af_mat_12 = P.*[-(rho_f1/tau1), -(rho_f2/tau1), -(rho_f1/tau1), -(rho_f2/tau1);
                  -(rho_f1/tau1), -(rho_f2/tau1), -(rho_f1/tau1), -(rho_f2/tau1);
                  -(rho_f1/tau1), -(rho_f2/tau1), -(rho_f1/tau1), -(rho_f2/tau1);
                  -(rho_f1/tau1), -(rho_f2/tau1), -(rho_f1/tau1), -(rho_f2/tau1)];
A_Af_mat_13 = (1/tau1)*eye(4);
A_Af_mat_21 = -kappa1*eye(4);

A_Af_mat_22 = [1-P(1,1)*beta1*rho_f1, -beta1*P(1,2)*rho_f2, -beta1*P(1,3)*rho_f1, -beta1*P(1,4)*rho_f2;
               -beta1*P(2,1)*rho_f1, 1-beta1*P(2,2)*rho_f2, -beta1*P(2,3)*rho_f1, -beta1*P(2,4)*rho_f2;
               -beta1*P(3,1)*rho_f1, -beta1*P(3,2)*rho_f2, 1-beta1*P(3,3)*rho_f1, -beta1*P(3,4)*rho_f2;
               -beta1*P(4,1)*rho_f1, -beta1*P(4,2)*rho_f2, -beta1*P(4,3)*rho_f1, 1-beta1*P(4,4)*rho_f2];
A_Af_mat_23 = zeros(4,4);
A_Af_mat_31 = [-gamma1*eye(2) zeros(2,2); zeros(2,2) -gamma2*eye(2)];
A_Af_mat_32 = [-alpha1*eye(2) zeros(2,2); zeros(2,2) -alpha2*eye(2)];
A_Af_mat_33 = eye(4);

A_Af_mat = [A_Af_mat_11,A_Af_mat_12,A_Af_mat_13;
            A_Af_mat_21,A_Af_mat_22,A_Af_mat_23;
            A_Af_mat_31,A_Af_mat_32,A_Af_mat_33];
 
        
solvec_f = inv(A_Af_mat) * [0;0;0;0;-kappa1;-kappa1;-kappa1;-kappa1;0;0;0;0];

A_Ai_mat = [1-p11*rho_i, -(1-p11)*rho_i, -(p11*rho_i)/tau1, -((1-p11)*rho_i)/tau1, 1/tau1, 0;
            -(1-p22)*rho_i, 1-p22*rho_i, -((1-p22)*rho_i)/tau1, -(p22*rho_i)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_i, -beta1*(1-p11)*rho_i, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_i, 1-beta1*p22*rho_i, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];        
solvec_i = inv(A_Ai_mat) * [0;0;0;0;1;1];


disp('')
disp('Solution: Coefficients on Technology Shock ')
disp(['output     : y_{1}^{a} (E1) = ' num2str(solvec_a(1))]); y1a = solvec_a(1);
disp(['output     : y_{2}^{a} (E2) = ' num2str(solvec_a(2))]); y2a = solvec_a(2);
disp(['inflation  : pi_{1}^{a} (B1) = ' num2str(solvec_a(3))]); pi1a = solvec_a(3);
disp(['inflation  : pi_{2}^{a} (B2) = ' num2str(solvec_a(4))]); pi2a = solvec_a(4);
disp(['interest rate : i_{1}^{a} (H1) = ' num2str(solvec_a(5))]); i1a = solvec_a(5);
disp(['interest rate : i_{2}^{a} (H2) = ' num2str(solvec_a(6))]); i2a = solvec_a(6);
disp('Solution: Coefficients on Markup Shock ')
disp(['output     : y_{1}^{f} (D1) = ' num2str(solvec_f(1))]); y1f = solvec_f(1);
disp(['output     : y_{2}^{f} (D2) = ' num2str(solvec_f(2))]); y2f = solvec_f(2);
disp(['output     : y_{3}^{f} (D3) = ' num2str(solvec_f(3))]); y3f = solvec_f(3);
disp(['output     : y_{4}^{f} (D4) = ' num2str(solvec_f(4))]); y4f = solvec_f(4);
disp(['inflation  : pi_{1}^{f} (A1) = ' num2str(solvec_f(5))]); pi1f = solvec_f(5);
disp(['inflation  : pi_{2}^{f} (A2) = ' num2str(solvec_f(6))]); pi2f = solvec_f(6);
disp(['inflation  : pi_{3}^{f} (A3) = ' num2str(solvec_f(7))]); pi3f = solvec_f(7);
disp(['inflation  : pi_{4}^{f} (A4) = ' num2str(solvec_f(8))]); pi4f = solvec_f(8);
disp(['interest rate : i_{1}^{f} (G1) = ' num2str(solvec_f(9))]); i1f = solvec_f(9);
disp(['interest rate : i_{2}^{f} (G2) = ' num2str(solvec_f(10))]); i2f = solvec_f(10);
disp(['interest rate : i_{1}^{f} (G3) = ' num2str(solvec_f(11))]); i3f = solvec_f(11);
disp(['interest rate : i_{2}^{f} (G4) = ' num2str(solvec_f(12))]); i4f = solvec_f(12);
disp('Solution: Coefficients on Monetary Policy Shock ')
disp(['output : y_{1}^{i} (F1) = ' num2str(solvec_i(1))]); y1i = solvec_i(1);
disp(['output : y_{2}^{i} (F2) = ' num2str(solvec_i(2))]); y2i = solvec_i(2);
disp(['inflation    : pi_{1}^{i} (C1) = ' num2str(solvec_i(3))]); pi1i = solvec_i(3);
disp(['inflation    : pi_{2}^{i} (C2) = ' num2str(solvec_i(4))]); pi2i = solvec_i(4);
disp(['interest rate : i_{1}^{i} (I1) = ' num2str(solvec_i(5))]); i1i=solvec_i(5);
disp(['interest rate : i_{2}^{i} (I2) = ' num2str(solvec_i(6))]); i2i=solvec_i(6);

